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1. Introduction 

The work is inspired by thermo-and photoacoustic imaging (see e.g. [21], [9l EH [HI [23] 
for some articles related to the subject), where the problem is the reconstruction of 
the absorption density from measurements of the pressure outside of the object. This 
is the Inverse Problems according to the forward problem, which maps the absorption 
density onto the pressure by solving the standard wave equation. Various reconstruction 
methods have been suggested in the literature for photoacoustic imaging, which can for 
instance be found in the excellent survey [10]. Recent efforts have been made to take 
into account attenuation [151 El E] and varying wave speed [6] . The standard model of 
attenuation (which is reviewed in Section [3]) is formulated in the frequency range and 
models the physical reality that high frequency components of waves are attenuated 
more rapidly over time. 

In this paper we review standard attenuation models, like power laws [TTj, [T9|, [20l 
[T8l[2l] and the thermo-viscous wave equation [8]. In this context, we discuss causality, 
which is the desired feature of attenuation models. The lack of causality of standard 
models in the parameter range relevant for photoacoustic imaging requires to investigate 
novel equations, which are derived in Section [3] and the following. 

We base the derivation of causal attenuation models on the mathematical concept 
of linear system theory, which can for instance be found in the book of Hormander [5j. 
In Section [H the abstract formulations are translated in equations which are formally 
similar to the wave equation. However, in general, the novel equations are integro- 
differential equations. An important issue is that the equations are formulated as 
inhomogeneous equations with homogeneous initial conditions, which is not standard for 
attenuated wave equations, where typically the equations are considered homogeneous 
with inhomogeneous initial conditions. For the standard wave equation these two 
concepts are equivalent, but only the one considered here, is mathematically sensible 
for the attenuation model. 

The approach leads to some novel causal attenuation models, in particular power 
law models (valid for a bounded frequency range), which are documented in the 
literature to be relevant for biological specimen (in the terminology used later on this 
means that 7 G (1, 2] - see [251 Chapter 7]) and also for instance also for castor oil, which 
satisfies a power law with index 7 = 1.66 [19]. These models are presented in Section [8] 
The rotationally symmetric examples, presented in Section [9l illustrate the unphysical 
behavior of some existing attenuation models. Aside from unmeaningly physical effects, 
the stable and convergent numerical implementation of attenuated, non-causal wave 
equations is an unconsidered problem since these equations lack the Courant-Friedrich- 
Levy (CFL) condition [^. The attenuation models considered here have a finite front 
wave speed and therefore can be implemented in a stable manner. Thus aside from 
physical considerations also from a point of view of stable numerical solution of wave 
equations questions of causality are most relevant. 

Concerning the presentation of the paper, the basic notation and mathematic results 
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are summarized in the appendix. 
2. Linear System Theory 

This section surveys linear system theory (see e.g. [13l [5]), which provides the 
hnk between hnear systems and convolution operators. This analysis is essential for 
the analysis. For notational convenience, when we speak about functions they are 
understood in the most wide meaning of the word, and can for instance be distributions. 
In the following, we give a characterization of causal functions and operators. 

Definition 2.1. (i) A function / := /(x, t) defined on the Euclidean space over time 
(i.e. in M^) is said to be causal if it satisfies 

/(x,t)=0 for t<0. 

(ii) In this paper A (with and without subscripts) denotes a real (that is, it is a mapping 
between sets of real functions on M^) and bounded operator. 

• ^ is translation invariant if for every function / and every linear 
transformation L := L{:K,t) := (x — Xo,t — to), with xq G and to ^ 

it holds that 

A{f o L) = {Af)L . 

Here o denotes the decomposition, i.e., (/ o L)(x, t) = /(L(x, t)) 

• ^ is called causal, if it maps causal functions to causal functions. 

• The operator A has a causal domain of influence if the function 

T(x) := sup{t : A6^^t{^, r) = for all t <t} for all x G 

is rotationally symmetric and the derivative with respect to the radial 
component T' satisfies 

< (T'ir))-' < Co < oo . (1) 
For convenience of notation we identify T(|x|) = T(x). 

The function T, presumably it exists, corresponds to the travel time of a 
wavefront initiated in at t = 0. ([1]) guarantees that the wavefront speed 
is finite. 

Remark 1. If the operator A models a physical process in a homogeneous and isotropic 
medium, then A is shift invariant and A6^^t is rotationally symmetric. 

If T exists, and in addition satisfies ([I]), then the property of a causal domain of 
influence guarantees that a wavefront can propagate with a speed of at most Cq. 

Now, we recall a fundamental mathematical theorem (see ^ Theorem 4.2.1]) of 
systems theory, which relates invariant operators with space-time convolutions. 
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Theorem 2.1. /5, Theorem Every linear (causal and) translation invariant 

operator A can be written as a space-time convolution operator with (causal) kernel 
G. That is, for arbitrary f from a suitable class of functions we have 

Af = G*^,tf. (2) 

In analogy to linear system theory we call the kernel G the Green function of A. 
According to Definition 12 . 1 1 the considered operators are real and therefore the according 
Green functions are real valued too. From the definition of the Green function it follows 
that 

G = A6^,t . 

In the following example we review the Green function and the convolution operator 
according to the wave equation. 

Example 2.1. We consider the standard wave equation in an isotropic medium with 
phase speed Cq G (0, oo): 



0. (4) 



together with initial conditions 

P\t<o = and —p 

'^^ i<0 

With source term / = 5x,t; the according solution Go of ([3]) and (jl]) is the Green function 

c 

47r Ixl 



5i t-ffl 

Go(x,t)= \ , . (5) 



Because of ([5]) Go is commonly denoted spherical wave originating from x = at time 
t = 0. 

In the space-frequency domain the Green function can be expressed by 

1 exp I luj— 
.F{Go}:=.F{Go}(x,^)= ^ 



V27r 47r |x 
It satisfies 



V.F{Go} 



lUJ 1 

.Co |x| 

and is the solution of the Helmholtz equation 



•^{Go}-sgn, (6) 



V'J' {Go} + — ^ {Go} = -^5x . (7) 

Cq V ZTT 

The operator 

-4.o/ := Go *x,t / 

is causal and maps a causal function / onto the solution of ([3]) and (jl]). 
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3. Attenuation 



In the chapter we investigated causality of attenuation models in a homogeneous, 
isotropic medium. In mathematical terms, it is common to describe attenuation by 
a multiplicative law in the frequency range: 

Definition 3.1. A real, bounded, linear, translation invariant operator A with causal 
domain of influence is called attenuation operator if there exists a complex function 
:= P^:{r,u) such that the associated Green function G := ^5x,t satisfies 



J^{G'}(x,cj) = exp(-/5,(|x|,a;))-JP{Go}(x,t<;) for all x G M^ G 



Here, JF is the Fourier transform (see Appendix). 

We rewrite ([8]) in the space-time domain by using 

1 



K := K{x,t) :-- 



/27r 



.F-i{exp(-/3.)}(|x|,t) . 



Therefore 



G'(x,t) = [K *t G'o](x,t) 



—K x,t- ^ 

Atx |x| \ Co 



(9) 



(10) 



Since in the context of this paper the operator A is real, the associated Green function 
is real- valued, and consequently [3^{r,u) has to be even with respect to u (cf. Property 
|v]in Appendix). 

Remark 2. In physical terms attenuation is a result of frequency dependent energy 
dissipation and therefore the ratio of the attenuated and un-attenuated wave amplitude 
must be smaller or equal to 1. That is 



exp (-5R(/3.)) 



T{G} 



< 1 . 



-^{Go} 

This implies that the attenuation coefficient (3^, satisfies 3?(/5=k) > 0. 

In the literature a special form of the attenuation coefficient is assumed: 

Definition 3.2 (Standard Form). The standard form of considered in the literature 
is (see e.g. [17] ) 

l3^{r,uj) = a^{uj)r for r > 0, u; G M. (11) 

The function a := ^{a^,) is called attenuation law. 

For the standard form P^, several properties for the attenuation operator are at 
hand. For instance the following results concerning travel time and causality. 

Theorem 3.1. Let A be an attenuation operator with (3^, of standard form. Then the 
travel time satisfies T(|x|) = |x| /c for some constant < c < Cq. 



Causality Analysis 3Tj^2s 6 

Proof. The definition of tlie travel time T in Definition 12.11 states that T(|x|) is the 
largest positive number such that for the Green function G = A6^^t 

G(x,t) = for t<r(|x|). 

This condition is equivalent to the condition that the function 

(x,t) G(x,t + T(|x|)) is causal (12) 

The operator A is causal and has causal domain of influence, which implies that T(0) = 
and {T'{r)y^ < cq. Consequently 



/•|x| 

T(|x|) = / r{r)dr > 
Jo 



(13) 



|x| 

Co ' 

r(r) := T(r) — r/co denotes the largest number such that -ft'(x, ■ + rdx])) = 
^^JF {exp (—/?,,)} is causal. From ([8]) it follows that 

This and the Theorem of Supports (cf. [5]) imply that r(r) = 2r(r/2), and consequently 
T is linear in r and after all T is linear as well. □ 

In particular, from ( |TOl) and Theorem 13.11 it follows that A has a causal domain of 
influence if and only if is a causal function. 

Remark 3. In the literature (for instance in [21]) causality is aimed to be enforced by 
demanding that 

JF~^ {a*} is causal. (14) 

This is equivalent to that the Kramers- Kronig relations for the m-th derivative ai"^^ of 
some function a* are satisfied, i.e., there exists a non-negative integer m, such that 

S>(a^))=7^{3fJ(a^))}=7^{a(™)} and := gfJ(a^)) = {53(a^))} , (15) 

where ?^ {•} is the Hilbert Transform (see Appendix). 

([T3j) follows already from the causality of K: From the definition of K it follows 
that ^ 

iV-f^l = —= \J^^^ {a* ■ exp {—a^ |x|)}| . 

Using some sequence {x^} with x„ 7^ and x„ — > shows that 

lim|Vir| (x„,t) = (t) . 

V 27r 



Due to the causality of K the left hand side is zero for t < 0, and thus is 
causal. 

However, as we show in Example 13.11 below, causality of JF^^ {a*} does not imply 
causality of K. In other words, in general, from the causality of JF^^ it cannot be 
deduced that A has a causal domain of influence. As a consequence several attenuation 
models considered in the literature lack causality. 
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Example 3.1 (Frequency power laws [IZllSn])- For the frequency power law, 



(16) 



where 7, ao > and 7 ^ N. The Kramers-Kronig relation with differentiation index 
m = 1 is satisfied for the one-parametric family of complex extensions (as considered 
[22 [20]) 

(17) 



cos 



(17) 



-100) ' + aoiio . 



Indeed O Theorem 7.4.3] implies that for every polynomial p in —iuj with nonegative 
real exponents, {p} is causal. Hence if 7 > 1, then 



(lu) := a^{uj) — ao {—iuj) Oq G 



has the same real part as and {o^f^} is causal. As a consequence the attenuation 
law a together with the causality condition (IT^ does not uniquely determine the 
attenuation operator A (cf. Definition 13. II) . 

Let be defined as in f|T7|) . then according to O Theorem 7.4.3] K, as defined in 
([9]), is causal if and only if 7 G [0, 1). Consequently, for frequency power laws with 7 > 1 
the according operator A, defined in Definition 13.11 does not have a causal domain of 
influence. 



4. Equations for Attenuated Pressure Waves 

In this section we formulate a causal wave equation which takes into account attenuation 
and review the literature (cf. [HI [201 [TH [IH [15]). 

Let A denote a translation invariant operator with causal domain of influence with 
travel time function T and cq as in Deflnition 12.11 The Green function G satisfles ([TOll 
and and therefore the according attenuation coefficient is given by 

A(x, u) = -log \^2y^{2^ |x| ^ |g (^x, . + I (^)| . (18) 

In the following we rewrite the term V^J-'{G} from which we derive the Helmholtz 
equation for JF{G}. Using f[TOl) . which states that G = K *t Gq = ^5x,t, and the 
product rule yields 

{G} = V^T {K} ■ T {Go} + 2V^ {K} ■ VT {Go} + T {K] ■ V^T {Go] ■ (19) 



To evaluate this expression, we calculate VJFjfC} and V^JF{K}. From ([9]) it follows 
that 

V^W = -/3:-^m-sgn, (20) 
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where denotes the derivative of (]^:{r,uj) with respect to r. This together with the 
formula fl62|) in the Appendix imphes that 

V^J^{K} = -V-{P:-J^{K}-sgn) 



- (V ■ sgn) ./3:.J^{K}- (sgn -VP'J-J^ {K} - (sgn ■ {K}) ■ (31 



(21) 

Inserting ([20]) and into and using the identity G = K *t Gq (cf. ([10])), shows 
that 



(22) 



- 2/5: ■ ^ m ■ (sgn ■ {Go}) + J" {K} ■ V^^ {Go} . 
Together with ([6]) and ([7]), the last identity simplifies to 

iu! 1 

Co |x| 



V2^{G} 



|X| 
,2 



J^{G}-2 



/3:-^{G} 



j-{G}-^-^{/0-5, 



(23) 



Cn 



/27r 



Since JF {i^} (x, u;) ■ (5x = ^ {K} (0, cj) ■ 5x, we obtain from ([23]) the Helmholtz equation 



V2^{G} 



/?: + 



Co 



^{G} 



/?:'-^{G}--=^{K}(0,cu)-5x 



^27r 



-/?:'■ ^ {G} - — exp (-/?, (0, u))-6^. 
zir 



(24) 



To reformulate ([24]) in space-time coordinates, we introduce two convolution operators: 



DJ := ir, / and D'J := /, 
where the kernels and i^'^ are given by 

1 

-fC* := -ft'*(x, t) := /C^,(|x| ,t) and K^{r,t) := - 
and 

:= i^':(x,t) := K{\^\,t) and if:(r,t) = -^T'' {/?:'} (r,t). 



/27r 

Using these operators and applying the inverse Fourier transform to ( [2^ gives 



(25) 



(26) 



(27) 




G = -D:G - K{0,t)6^. 
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For a general source term /, patt := -^f = G *x_t / solves the equation 



V72 1 



-Asf 



(29) 



where As denotes the space-time convolution operator with kernel 

Ks := i^.(x, t) := -BG + D',G + K{0, t) ■ 6^ 

where 

Co at 

Equation ( l29i) is the pressure wave equation that obeys attenuation with attenuation 
coefficient ffTSl). 



(30) 



(31) 



Remark 4. In this remark we consider again the standard model, as in Definition 13.21 
For this case the wave equation fl29|) can be casted in a form that resembles the standard 
attenuation wave equation (cf. Example 14.11) . Since K is causal, it follows that K^, is 
causal too (the argumentation is analogous to Remark [3]) and therefore the operator 
is well-defined for all causal functions. Moreover, since K'^ = 0, it follows that D'^ = 0. 
Using that depends only on t it follows that 

(D^G) / = [K* *t G] f = K,*t [G *x,i /] = D,{G *x,t /)• 

Convolving each term in (!28l) with a function /, using the previous identity and that 
= 0, it follows that 



l_d_ 

Co dt 



Patt — ~/ • 



(32) 



In the following we review some wave equations obeying attenuation, which are 
frequently considered in the literature: 

Example 4.1. • For 7 > and 7 ^ N, denote by D] be the Riemann-Liouville 
fractional derivative (see OUn])) with respect to time. It is defined in the Fourier 
domain by 

J^{D^f} = i-iujrj^{f} (33) 

and satisfies 



Dpf = D2D2f 



and 



d 

dt^'^ 



d 



Now, we consider the attenuation coefficient 

/3*(r, cij) := aoi^—iuyr with ao := ao/ cos(77r/2) 
which satisfies the attenuation law 

^{P*){r,u) = a{u!)r and a^u) = ao\uj\'^ 



(34) 



(35) 
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(cf. Example 13.11 and 
with kernel K^, defined by 



). Let denote the time-convolution operator 
and fl35l). Then form fl33l) and = 



T-^ {ao(-i^)^} /V27r it follows that = a^D^. In (mill] (see also [IHlEn!) the 
following equation for the pressure function patt of attenuated waves is investigated: 



(36) 




which is equivalent to equation (!32|) with operator D^, = ckoDj. Let A denote the 
solution operator of (l36ll . then from [5^ Theorem 7.4.3] it follows that A has a causal 
domain of influence only for 7 G [0, 1). 

Let 7 > 0,7 ^ N. Neglecting in fl36l) the operator a^Dp (which one finds after 
expanding the decomposition operator) one finds Szabo's equation [T9] 



V Patt - TfTil^^'att 



Co 



(37) 



This equation is equivalent to equation fl32|) if we define the kernel of by fl26l) 
with 



/3*(r, tu) := Q;o(-ia;)'''r and a^{ijj) 



Co Co 



{-lujf + 2aoCo(-icu)^+i 



(38) 

Again, if A denotes the solution operator of (|57|) . then [51 Theorem 7.4.3] implies 
that A has a causal domain of influence only for 7 G [0, 1). 

In the literature, the standard attenuation models ( l36l) and ( l37l) are considered 
as homogeneous Cauchy problems with inhomogeneous initial conditions. In contrast, 
in our setting, we consider inhomogeneous Cauchy problems with homogeneous initial 
conditions. In the following section we show that the two concepts can be equivalent. 
However, in general, only the concept suggested here leads to a rigorous framework, in 
which we can define solution operators for attenuated wave equations. 

For the readers convenience, we summarize some important notation and facts in 
the following table. Note the difference between K, K^, and K'^, respectively, with 
respect to the involved exponential function. 



Kernel 


General 


Standard Form 


Convolution Operator 


K ^ 


^^-i{exp(-/5.)} 


^-^"Uexp(-«* |x|)} 




(EE]) 








K m 







d: m 



5. The homogeneous Cauchy problem with memory 

We consider the standard attenuation model /5^,(r, cj) = a^:{uj)r. Let A denote a 
translation invariant operator with causal domain of influence and and let the operator 
be as defined as in fl^ . 
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In this section we investigate under which conditions the inhomogeneous wave 
equation (|32|) with homogeneous initial conditions (jl]) (where p is replaced by patt) 
and the homogeneous equation 



V'gatt 



]_d_ 

Co dt 



gatt = 



(39) 



with the inhomogeneous initial condition 

Qutt = qo for t <0 and 



d_ 

di 



t=0+ 



d_ 

di 



(40) 



t=o- 



are equivalent. That is, both equations have the same solution for t > 0. 



Theorem 5.1. Assume that [3^) . and [3^) . have unique solutions, respectively. 
Then qatt = Patt for t > if and only if qo and patt are related by the following conditions 



^p:= hm = lim w:Patt 
t^o- dt t^o+ dt 



and 



1 / d 
H ■ Bqo = -f + [ip ■ 6t + (f ■ -^dt j , 



(41) 



(42) 



with B is as in [3l\) and H is the Heaviside function. 

Proof. Assume that gatt = Patt for t > 0. Then, using that Patt = for t < 

implies that 

Patt = H ■ gatt and gatt = Patt + qo ■ 



In particular, property (HT!) holds. Moreover, (H3ll implies 



V^Patt = H ■ V^gatt and 



92 9^ d 

Q^P-tt = H ■ —gatt + ^ ■ <5t + <^ ■ ^5^. 



Since 



+ 

cldt^ 



Co dt 



(43) 



(44) 



(45) 



it follows from (HI, (|45]) and (|43]) that 



V'patt - Bp 



1 d^ 

-2^^ Patt + f = H 



c^dt^ 



2 R 1 92 ■ 

V gatt ~ ^Q'att oTTT^ Qatt 

Cq dt^ 

Bpsxt + H ■ Bq^tt - ^ f ^ ■ 5t + ■ ) + / • 



dt 



(46) 



Using the definitions of gatt and patt, ( 14611 simplifies to 

d_ 
di 



-Bp,tt + H ■ Sgatt -^(^■6t + v-^Jt]+f = 0. 
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Since is a causal operator and Patt is a causal function, we have Spatt = H ■ Spatt- 
This together with P3l) implies that — iSpatt + -f^ ■ -S^att = -f^ ■ ^Qo- Hence 

H-Bqo-^^(^iJ-St + ^- ^^S?j + / = 

and thus (l42l) holds. This proves the first direction of the theorem. 



=: To prove the opposite direction let 

Patt := H ■ gatt such that 



gatt = Patt + go • 



We prove that Patt = Patt holds for t > 0. Similarly as in part a) of the proof it 
follows that 



V Patt ^ ^'Patt 9 ol9 Patt 



2 1 92 ■ 

V (Jatt ~ S^att 2 "0I9 5'att 



1 ( d 



holds. Since gatt solve problem (jSlD, (HO]) and condition 
identity simplifies to 



is satisfied, the last 



V Patt ^'Patt yWTo Patt 



Hence we have shown that patt solves problem ( 132|) . (jll) and since this problem has 
the unique solution Patt, it follows Patt = Patt- In summary we have shown that 



Patt — Patt — ^att 

which proves the assertion. 



for 



t > 



□ 



Remark 5. In the absence of attenuation the operator B is the zero operator and 
condition fj42l) reduces to 

In this case the solutions of 



dt 



V^p- 


1 92 „ _ f 

c^dt^P - J ) 




1 9\ - n 


P\t<o 


= = 


glt=o 





are identical for t > 0. 



Causality Analysis 3T£]X2s 

6. The thermo-viscous wave equation 



13 



In this section we show that the thermo-viscous wave equation (see e.g. [8]) is not causal 
(see Theorem 16.1 1 below) . The formalism introduced here will enable us to derive a causal 
variant of the thermo-viscous equation which satisfies the same attenuation law. 

The thermo-viscous wave equation models propagation of pressure waves in viscous 
media and reads as follows 



-F . 



(47) 



Here tq and cq denotes the relaxation time and the thermodynamic speed, respectively 
and F models sources. 

In the following we transform the thermo-viscous wave equation into the form fl28|l . 
which enables us to deduce that the thermo-viscous equation is not causal. For these 
purpose we consider the attenuation coefficient 



/3^,(r, Co") = a^{uj)r 



with 



la; 



Co Co VI - Toi^ 
and the time convolution operators T^^^ and L^^^ with kernels 

Kj,,/2 := -^T-^ {(1 - irocu)-^/^} and 
V 27r 



(4^ 



V ZTT 



respectively. Since K^, satisfies (1261) it can be rewritten in the following form 



1 ^ 1 \ lUJ 



^2tc I Co Co a/I - i^^^o 
Therefore the according convolution operator D^: is given by 

1 d ^ d 



Co dt * Co dt '^^''^ ' 



= — 



dt 



(49) 



(50) 



Co dt Co ' 

In the following we summarize some properties of the operators T^/^, L^/^ and D^, and 
the associated kernels. 

Lemma 6.1. The kernel functions Krpi/2, and the operators L^^^, T^^"^ , 

respectively, satisfy: 

(i) For To = 0, = = I and = /sT^i/a = 6t. 

(n) 



KT.,2{t) 



^H{t) exp(-t/ro) 



r(l/2)ro/V/2 



(51) 



(Hi) L}/'^ is the inverse ofT^^'^, 
(iv) Let L := [L^l^f and T := [T^'^f , th 

T T 9 



en 



and 



T = L-\ 



(52) 
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(v) D[ = and for tq = also = 0. 



Co at 



Proof, (i) The first item is a trivial consequence of properties of the Fourier transform 

:f-^{-}. 

(ii) With the substitution s = —iujtq we derive the relation with the inverse Laplace 
transformation {■} (for a definition and some basic properties see the appendix 
of this paper) 

T-^ {(1 - i^ro)-^/'} (t) = : ^ / exp (st/ro) ■ (1 + sY^I'^ds 

IToV 27r J-ioo i^'i^ 



Using the properties (164|) and (!65|) of the inverse Laplace transformation the 
assertion follows. 

(iii) From 

V 2 7r 



(54) 



it follows that for each function / 

^1/2 ^ ^ ^^^^^ ^^^^^ ^^f = S,*,f = f 
^1/2 2^1/2 ^ ^ ^^^^^ ^^^^^ ^^f = S,*,f = f. 

(iv) Since 

1 _ 9 

Kli/2 *t Kli/2 = —=T~^ {1 - icJTo} = (5f - ro^(5t, 
V ivr Ct 

it follows that 

The assertion T = is then a consequence of the previous item. 

(v) Since K^: does not depend on |x| and K'^ is the kernel of D[, it follows that K'^ = 0, 
i.e. D[ = 0. The second statement is a direct consequence of Item [i] which states 
that = / for tq = 0. 

(vi) Follows from (|50l) . 

□ 

The thermo-viscous wave equation (H71) can be put in formal relation to the wave 
equation fl28l) by identifying an appropriate operator as in fl50|) : 
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Utilizing Item|vI]of Lemma ED in equation fl32p and taking into account that D'^ = 
(cf. Item 13 of Lemma [6. ip shows that the solution of the thermo- viscous wave equation 
(HTl) with F := Lf satisfies 



VVatt - 



]_d_ 

Co dt 



1 2 



Patt 



2 1 



Conversely, the solution of equation ( 1321) with defined as in (1501) satisfies the thermo- 
viscous wave equation fj47|) with F = L f. 

Theorem 6.1. Let a* be defined as in ( f^^ . T/ien 3?(a*) anc? 53(q;=k) satisfy the Kramers- 
Kronig relation, hut the solution operator A of the thermo -viscous wave equation does 
not have a causal domain of influence. 

Proof. Since defined as in fHOl) is causal, it follows that '^{a^) and 55(q;*) satisfy 
the Kramers-Kronig relation. From [5l Theorem 7.4.3] it follows that the kernel 
K := -^^T^^ {exp (—a* |x|)} is not causal and as a consequence the according solution 
operator of the t her mo-viscous wave equation does not have a causal domain of 
influence. □ 

Remark 6. From (j48|l it follows that the attenuation law a = ^{a^,) approximates for 
small frequencies the frequency power law with 7 = 2. 



7. A causal thermo-viscous wave equation 

Below we discus a causal variant of the thermo-viscous wave equation. 

Let «! > 0. Theorem 16.11 below shows that the attenuation operator with 
attenuation coefficient of standard form P^{r,u!) = a'^{u)r and 



aiiuj 



CoVl - To'lUJ 

has a causal domain of influence. The operator and its kernel i^^, read as follows 



(55) 



Co 



dt 



and 



a\ d 
Co dt 



K. 



(56) 



Note that D'^ = 0, since K^, does not depend on |x|. 

For «! = 0, = and thus equation ( l28l) with operator D^, defined by ( l56l) is the 
standard wave equation (without attenuation). Since 

and L = 



2 ^ 
dt^ 



281) can be rewritten as 



To 



d_ 

di 



1^ 



Patt 



d_ 
di 



(57) 
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Theorem 7.1. Let and a% he defined as in and / f53]) . respectively. Then 



3fJ(Q!^) and Q{al) satisfy the Kramers- Kronig relation and 3fJ(a*) = 3?(a^). The solution 
operator A of equation [5l\ ) has a causal domain of influence. 

Proof. Since K^, defined as in fl56l) is causal, it follows that ^{a^) and 3(a;^) satisfy 
the Kramers-Kronig relation. Comparison of defined as in ( l48l) and defined as 
in (155!) shows that 3?(a*) = 3?(q;^). From [5l Theorem 7.4.3] it follows that the kernel 
K := -^^'^ {exp (— |x|)} is causal and as a consequence the solution operator of 
equation fl571) has a causal domain of infiuence. □ 



Remark 7. In ultrasound imaging soft tissue is often modeled as a viscous fiuid and 
therefore fl57j) is a potential model, on which thermoacoustic tomography can be based 
on. Moreover, the attenuation of tissue is frequently modeled as a power frequency law 
with 7 G (1,2). 

8. Causal Wave Equations satisfying Frequency Power Laws for small 
frequencies with 7 G (1,2] 

In Example 13.11 we have shown that the frequency power law does not yield to a 
causal wave equation when 7 > 1. In this section we derive causal wave equations for 
attenuation laws which approximate frequency power laws for small frequencies with 
exponent 7 G (1,2], where for 7 = 2 we get the causal variant of the thermo- viscous 
wave equation (1571) . 

Here we follow the notation of the previous section and introduce the following 

1/2 1 /2 

families of operators: For constants 7 G (1,2], tq > and ai > let and L-y 
denote time convolution operators with kernels: 

K^y. := 4=-^ I (1 + i-iooToy-T'^'] , 

^= vb•^{(^ + (-^"^°)'")'''}• 
We set := (r^^^y and := (l^^)^. We emphasize that T^^^ = T^/^ where T^/^ 

is the operator in the thermo-viscous case. The operators T^^'^ and satisfy similar 
properties as the operators T^/^ and L^/^ in the thermo-viscous case: 
The following lemma is proven analogously as Lemma 16.11 

Lemma 8.1. • For t = we have 

T^/^ = L\l^ = I and = K^y^ = 5t . 

• is the inverse ofT^^"^. 

• Let -D7~^ ^6 fractional derivative of order 7 — 1, as defined as in [33^) . then 
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In analogy to Section E] we consider now the standard attenuation coefficient with 

a*H = , . . =^- (58) 

Here ai, Tq and Cq are positive constants that are medium specific. The operator D^, 
and its kernel fC^, are given by 



:= -n"^^ and = ^^K^.;. • 

Co ' at Co at -'t 



(59) 



Moroever, the kernel K, defined by ([9]), reads as follows 

(t) (60) 



i^r(x, t) = ^— J^ I exp ( , Q;ii^|x| 

\ \cov/l + (-iroCc;)^-\ 



For uj small we have 

Q;(u;j ~ ao \touj\ ' with Oq = . 

2roCo 

The wave equation fl28|) with as in fl59|) reads as follows 

(/ + rj-^ A"-^) V Vtt - [ai / + Lf ] ' = -F. (61) 



c, 



In particular, for 7 = 2 we recover the causal variant of the thermo- viscous wave equation 

(E 



Theorem 8.1. The solution operator of equation l[61\) has a causal domain of influence. 

Proof. From [Sj Theorem 7.4.3] it follows that K from (!60|) is causal and thus the solution 
operator of (16T|) has a causal domain of infiuence. □ 

9. Examples 

In this section we present some calculations, highlighting the effects of non-causality. In 
all examples is of standard form ( ITTl) and the solution operator A determined by 
has the Green function 

1 r exp (-/?*) ■ exp ( -i ■ (t - g) ) 
V2^Jm 47r|x| 

We recall that the operator A has a causal domain of infiuence if and only if 
J-""^ {exp (— /5=k)} is a causal function. In other words, non-causality can be observed 
if 

{exp (-/9*)} (t) ^ for some t < . 
All numerical simulations were performed in MATLAB with the fit-subroutine. 
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Figure 1. Simulation of T ^ {exp (— /3»(|x| ,uj))} for the frequency power law with 
(7,ao) e {(0.5,0.1581), (1.5,0.0316), (2.7,0.0071), (3.3,0.0027)}, cq = 1 and |x| = i. 
In the first example 7 < 1 and thus the function is causal. For all other cases it is non 
causal as predicted by the theory. 

Frequency power law: Let a = ao lul^ with some 7 > 0. The extension 
by the Kramers-Kroenig relation is given by (JTTIl . Fig. [1] shows simulations of 
JF~^ {exp (— /?*)}, which illustrates that causality only holds for 7 G [0, 1). 

Szabos's model: Here a^:{u!) is as in (!38|) . In Fig. [2] we show simulations of 
{exp (— /5*)}. The numerical result confirm the mathematical considerations 
that causality only holds for 7 G [0, 1). 

Thermo- viscous wave equation: There a^, is as in fHSl) . The left pictures in Fig. [3] 
shows a simulation of {exp (—/?*)} for the thermo- viscous wave equation (H7|) . 
Note that according to (HSj) and (1551) the attenuation laws of the thermo-viscous 
wave equation and the causal variant (|57|) differ just by a multiplicative constant 
«i. A simulation of JF"^ {exp (— with ai = 1 for the causal variant (!57l) of the 
thermo-viscous wave equation is shown in the right pictures of Fig. [3] 

10. Conclusions 

In this paper we introduced the concept of an operator with causal domain of influence 
which guarantees a finite wave front speed. As a consequence these models allow for 
a stable numerical implementation and thus are suitable for photoacoustic imaging, 
where inversion techniques are required. Based on this concept, we showed that an 
attenuated wave described by such an operator satisfies the standard causality condition 
known as the Kramers- Kronig relation. However theses relations are not sufficient to 
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guarantee that an attenuated wave has a finite wave front speed. This is a common 
misunderstanding in causahty theory. 

We also showed that attenuated waves satisfying the frequency power law and the 
Kramers-Kronig relation have finite wave front speed only if 7 G (0, 1). An example of an 
equation where waves can propagate with infinite wave front speed is the thermo- viscous 
wave equation. Because of the lack of causality of standard models in the parameter 
range relevant for photoacoustic imaging, we developed novel equations that satisfy our 
causality requirement and the desired attenuation properties. 

For our causality analysis all equations were formulated as inhomogeneous equations 
with homogeneous initial conditions, but we showed that if certain conditions are 
satisfied, then the attenuation problem can be formulated as a Cauchy problem with 
memory. 

11. Appendix: Nomenclature and elementary facts 

Real and Complex Numbers: C denotes the space of complex numbers, R the space of 
reals. For a complex number c = a + ib a = 9fJ(c), b = Q{c) denotes the real and 
imaginary parts, respectively. 

Differential Operators: V denotes the gradient. V- denotes divergence, and denotes 
the Laplacian. 

Product: When we write ■ between two functions, then it means a pointwise product, 
it can be a scaler product or if the functions are vector valued an inner product. 
The product between a function and a number is not explicitly stated. 




Figure 2. Simulation of ^ {exp (— /3*(|x| , w))} for Szabo's frequency law with 
(7,ao) e {(0.5,0.1581), (1.5,0.0316), (2.7,0.0071), (3.3,0.0027)}, cq = 1 and |x| = i. 
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Figure 3. Left: T^^ {exp (— /3, (|x| , lS))] defined by the thcrrno- viscous wave equation 
(|Tf)) with To = 10"^, Co = 1 and fixed |x| = \. Right: Causal variant ([57]) of the 
thermo-viscous wave equation with ai = 1, tq = 10^^, cq = 1 and fixed |x| = j. 



Decomposition: The decomposition of operators A and B is written as AB. 
Special functions: The signum function is defined by 



sgn := sgn(x) := - — - 



In it satisfies 



2 

V ■ sgn = — . (62) 



X 



The Heaviside function 

H := H{t) 

satisfies 



for t < 

1 for t > 



H := ^(1 + sgn) . 

The Deha-distribution is the derivative of the Heaviside function at and is denoted 
hj St := St{t). In our terminology St denotes a one-dimensional distribution. The 
three dimensional Delta-distribution S^ is the product of three one-dimensional 
distribution Sx^, i = 1,2,3. Moreover, 

S^,t ■= 5x,t(x, t) = S^-St, (63) 

is a four dimensional distribution in space and time. 

Properties related to functions: supp((7) denote the support of the function g, that is 
the closure of the set of points, where g does not vanish. 

Derivative with respect to radial components: We use the notation 

r := r(x) = |x| , 

and denote the derivative of a function /, which is only dependent on the radial 
component |x|, with respect to r (i.e., with respect to |x|) by ■'. 
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Let f3 = f3{r), then it is also identified with the function (] = /3(|x|) and therefore 

V/? = . 
|x| 

Convolutions: Three different types of convolutions are considered: *t and denote 
convolutions with respect to time and frequency, respectively. Let /, /, g and g be 
functions defined on the real line with complex values. Then 

f*t9--= f fit - t')g{t')dt\ f *^g:= [ f{uj - u')g{uj')du;' . 



*x,t denotes space-time convolution and is defined as follows: Let /, g be functions 
defined on the Euclidean space with complex values, then 

f*^,tg-= / / f{-i^-:i^,t-t')g{-x!,t')diLdt . 

Fourier transform: For more background we refer to [22| [T3| [5]. All along this 
paper T {•} is the Fourier Transformation with respect to t, and the inverse Fourier 
transform JF^^ {■} is with respect io uj. In this paper we use the following definition 
of the Fourier transform T {■} and its inverse JF^^ {■} 

^{/}H = ^ / exp(ia;t)/(i)dt, |/| (t) = ^ / exp (-ia;t)/Hrfu; . 

V 2,71 Jm. ^ V 2,71 Jr 

The Fourier transform and its inverse have the following properties: 



(ii) 



^{|/|^H = (-M-^{/}M- 



^{/•^?} = ^^{/}*^^{(7} and 

V 27V 



J" {f} ■ J" {g} = {f *t g} 



T-^ {f-9\ = {/} ** and 



(iii) For a e M 

^ {fit - a)} {u) = exp {-lau) ■ ^ {fit)} (cu) 

(iv) The Delta-distribution at a G M satisfies 

Stit — a) = {exp(iaco')} (t) . 

(v) Let / be real and even, odd respectively, then JF {/} is real and even, imaginary 
and odd, respectively. 
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The Hilbert transform for L^— functions is defined by 

n{f}it) = --f l^ds, 

where j-j^f{s)ds denotes the Cauchy principal value of J^f{s)ds. 

A more general definition of the Hilbert transform can be found in [1] . The Hilbert 

transform satisfies 

• ^{^{/}}M = -2^{sgn/}M, 

. n{n{f}} = -f. 

From the first of these properties the Kramers-Kronig relation can be formally 
derived as follows. Since f{t) is a causal function if and only ii f = H ■ f and 
H = (l+sgn)/2, it follows that J^{/} = [J^ {J^ {/}}]/2, which is equivalent 

to^{/} = z7^{^{/}}, i.e. 

= and {f}) = {/}}). 

The inverse Laplace transform of / is defined by 



for t < 0, 

tj:-::^^vi^t)mds, for t>o, 



where 7 is appropriately chosen. 

The inverse Laplace transform satisfies (see e.g. [1]) that 

£-1 {h{s - a)} (t) = exp (at)£-^ {h{s)} (t) for all a,teR (64) 

and 

{s-^} (t) = (r > 0) . (65) 
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